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We study Mott transition in the two-dimensional Hubbard model on an anisotropic triangular 
lattice. We use the Lanczos exact diagonalization of finite-size clusters up to eighteen sites, and 
calculate Drude weight, charge gap, double occupancy and spin structure factor. We average 
these physical quantities over twisted boundary conditions in order to reduce finite-size effects. 
We find a signature of the Mott transition in the dependence of the Drude weight and/or 
charge gap on the system size. We also examine the possibility of antiferromagnetic order 
from the spin structure factor. Combining these information, we propose a ground-state phase 
diagram which has a nonmagnetic insulating phase between a metallic phase and an insulating 
phase with antiferromagnetic order. Finally, we compare our results with those reported in the 
previous theoretical studies, and discuss the possibility of an unconventional insulating state. 

KEYWORDS: Mott transition, geometrical frustration, Hubbard model, Lanczos exact diagonalization, 
Drude weight, averaging over twisted boundary conditions 



1. Introduction 

Mott transition is a metal-insulator phase transition 
driven by strong Coulomb repulsion 1 and has been a cen- 
tral theme in physics of strongly correlated electrons. 2 A 
Mott insulator is a featureless insulating state with an 
odd number of electrons per unit cell and can be distin- 
guished from a band insulator which has an even number 
of electrons per unit cell. In most cases, however, would- 
be Mott insulators show some (translational) symmetry 
breaking at low temperatures, resulting in a change in 
the size of the unit cell. For example, the Hubbard model 
and its generalizations defined on a bipartite lattice show 
antiferromagnetic order in their insulating phase at half 
filling (one electron per site). In the narrow definition 
of a Mott insulator, such insulating states are not Mott 
insulators but rather band insulators, as the unit cell is 
doubled and contains two electrons. It is hard to find a 
true Mott insulator without any symmetry breaking in 
nature. 

A promising route to a Mott insulator is to sup- 
press symmetry breaking by introducing strong geomet- 
rical frustration into lattice structure. A series of lay- 
ered organic conductors, ^-(ET^X, have recently at- 
tracted much attention as materials that may realize 
this route. 3 ' 4 The electronic properties of these com- 
pounds are determined by conducting carriers moving 
on an anisotropic triangular lattice, which is geometri- 
cally frustrated; the degree of frustration is controlled 
by the choice of the anion X. Furthermore the strength 
of electron correlation can be easily tuned by applying 
pressure. 

The k-(ET)2X compounds show various intriguing 
properties. For example, k-(ET)2Cu[N(CN)2]C1 shows a 
first-order Mott phase transition without any symmetry 
breaking at finite temperatures under pressure. In the 
vicinity of the critical end point of the first-order tran- 



sition line, the conductivity was found to show an un- 
conventional critical scaling behavior, 5 which is different 
from the Ising universality class expected to describe the 
Mott criticality in higher (than two) dimensions. 6-8 An- 
other compound k-(ET)2Cu2(CN)3, whose lattice struc- 
ture is close to the isotropic triangular lattice, has an in- 
sulating phase at ambient pressure showing no magnetic 
order down to very low temperature, T = 32mK. 9 ~ n 
Some suspect that this insulating phase is a Mott in- 
sulator with a spin-liquid ground state. When pressure 
is applied, the Mott insulating state undergoes a phase 
transition to a superconducting state. 

A minimal theoretical model that contains essential in- 
gredients for describing the electronic properties of these 
organic compounds is a single-band Hubbard model on 
an anisotropic triangular lattice at half filling. 12 It is 
thus important to understand the Mott transition and 
the resulting Mott insulating state in this and related 
models. Indeed there have been several theoretical stud- 
ies along this line. In particular, the dynamical mean 
field theory (DMFT) has been successfully applied to 
the study of the Mott transition in the Hubbard model 
at finite temperatures. Among other things it has pre- 
dicted that the Mott transition at the critical end point 
should belong to the Ising universality class. 6 However, 
the DMFT is an approach which is justified only in large 
spatial dimensions, and it cannot be directly applied to 
the two-dimensional case. Some extensions of the DMFT 
to two dimensions, such as a cellular DMFT, have also 
been applied to this model and led to results suggest- 
ing the existence of a first- or second-order transition 
without symmetry breaking at finite temperatures. 13-15 
However, it is not clear whether spatial correlations are 
properly taken into account in this type of approach. 
Other theoretical approaches to finite-temperature Mott 
transitions include a phenomenological effective theory 16 
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and a mean- field theory, 17 which pointed out an impor- 
tant role played by a marginal quantum critical point in 
the unconventional critical behavior near the critical end 
point in k-(ET) 2 Cu[N(CN) 2 ]C1. 

The possibilities of the ground state being a non- 
magnetic insulating state and a superconducting state 
have been studied for the Hubbard model on an 
anisotropic triangular lattice using many different ap- 
proaches such as the fluctuation exchange approxima- 
tion, 18 ' 19 a U(l) gauge theory, 20 variational Monte Carlo 
(VMC) studies 21 ' 22 and the cellular DMFT. 23 Related 
models were also discussed, for example, in the varia- 
tional study of the Heisenberg spin model with ring ex- 
change coupling 24 and the resonating- valence-bond mean 
field theory for the Hubbard-Hciscnbcrg model. 25,26 
These studies, however, often rely on some uncontrolled 
approximations whose validity is uncertain when both 
the strong correlation and the geometrical frustration are 
important. 

For the half-filled Hubbard model with geometrical 
frustration, choices are quite limited of numerical ap- 
proaches that can both treat correlations in real space 
and work at zero temperature (or low temperatures). 
For example, the DMFT and its cluster extensions do 
not fully treat the correlations in the real space, as we 
have mentioned. Monte Carlo simulations are suscepti- 
ble to the infamous negative sign problem, and work only 
at high temperatures. At this stage, available numerical 
methods that satisfy the required conditions are the path 
integral renormalization group (PIRG) method 27 ' 28 and 
the exact diagonalization method. 

Imada and his coworkers have developed the PIRG 
method and applied it to study the Hubbard model on 
two-dimensional frustrated lattices. In the ground-state 
phase diagram of the half-filled Hubbard model on an 
anisotropic triangular lattice, they found a narrow re- 
gion of a non-magnetic insulating phase between a para- 
magnetic metal and an antiferromagnetic insulator, if 
the frustration is substantial. 29 In the non-magnetic in- 
sulating state, no obvious symmetry breaking pattern 
was observed, and energy spectrum has highly degen- 
erate low-lying excitations as in the Fermi liquid, which 
could be attributed to the presence of a spinon Fermi sur- 
face. 20 ' 24 The transition from the metallic phase to the 
non-magnetic insulator appears to be first order although 
the discontinuity becomes weaker for strong frustration. 
In order to confirm and better understand these inter- 
esting behaviors, a complementary study using another 
unbiased technique is needed. 

To this end, we use the exact diagonalization method 
to study the ground-state properties of the two- 
dimensional Hubbard model on the anisotropic triangu- 
lar lattice at half filling. The method treats correlation 
effects in an unbiased way with high precision, although 
the calculation is limited to small systems. We exam- 
ine dependence on boundary conditions (BCs) for finite- 
size clusters, and demonstrate that averaging over the 
twisted BCs is an efficient tool to reduce the finite-size 
effects. We calculate the Drude weight, which was not 
calculated in the previous works, to monitor the metal- 
insulator transition directly. We also compute the charge 
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Fig. 1. Anisotropic triangular lattice and the hopping integrals 
in our model (1). 



gap, the double occupancy and the spin structure fac- 
tor to elucidate the phase diagram and the nature of the 
transitions. The results are discussed in comparison with 
those in the previous works. 

This paper is organized as follows. In §2, we define the 
model and finite-size clusters to be studied numerically. 
The physical quantities signalling the Mott transition are 
also introduced. In §3, we discuss their dependences on 
BCs in the finite-size clusters, and demonstrate that av- 
eraging over the twisted BCs leads to substantial reduc- 
tion of the finite-size effects. In §4, we present numeri- 
cal results on the ground-state properties. The numerical 
results and their implications to the phase diagram are 
discussed in §5. We conclude with a brief summary in §6. 

2. Model and Method 

We study the ground-state properties of the Hubbard 
model on an anisotropic triangular lattice at half filling. 
Its Hamiltonian is given by 

H=-J2 kjixijclrCje + H.c.) + [/^n, T n a , (1) 

{ij)o i 

where Ci a (c\ a ) is the annihilation (creation) operator for 
electrons of spin a on the site i, and rii a — c\ a c i(y is the 
number operator. The model is defined on the anisotropic 
triangular lattice with the hopping integrals tij as shown 
in Fig. 1. We take t = 1 as the energy unit and fix the 
electron density at half filling, one electron per site on 
average. 

We consider finite-size clusters on the triangular lattice 
shown in Fig. 2. To impose a twisted BC for each cluster, 
we have introduced the phase factor, 

Xij =exp(i0-r ij ), (2) 

where is the vector connecting from the site i to the 
site j, and 4> = 4>\b\ + ^> 2 &2 in terms of the reciprocal 
vectors foi and b 2 , which satisfy the orthonormal rela- 
tions di ■ bj = Sij with the vectors a\ and a 2 defining 
the cluster (Fig. 2). The "flux" <p = (<j> u (f> 2 ) dictates the 
BC along the direction parallel to the vectors a\ and 
a 2 . For example, <f> = (0,0) represents periodic BCs in 
both directions, and <f> = (0, w) periodic BC in the a\ 
direction and anti-periodic BC in the a 2 direction, re- 
spectively. Any nonvanishing flux cf> leads to a twisted 
BC (-7T <<j>i<n,i = 1,2). 
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Fig. 2. Finite-size clusters used in our study. The number indi- 
cates the size of each cluster, N. There are two different geome- 
tries for the 12-site clusters. All the clusters are compatible with 
antifcrromagnetic Neel order while only the 12(a)- and 18-sitc 
clusters are compatible with 3-sublattice order. 



We apply the exact diagonalization method based on 
the Lanczos technique to finite-size clusters up to 18 
sites shown in Fig. 2. Since the model has translational 
symmetry, each energy eigenstate is classified according 
to its total momentum k. For a given </>, we numeri- 
cally obtain a ground-state wave function of a finite-size 
cluster for several values of U and t'. To save the CPU 
time, for 16-site and 18-site clusters, lowest-energy states 
are computed at some relevant momenta only, such as 
k = (0, 0), (ir, 7r), (it, 0), (0, 7r), and the momentum of the 
ground state at U = 0. These wave functions on finite- 
size clusters depend on 4>. As we discuss in the next sec- 
tion, the ^-dependence can be used for reducing finite- 
size effects. Here let us assume for the moment that <fi is 
fixed at a certain value. 

To study the metal-insulator transition at zero tem- 
perature, we compute the Drude weight and the charge 
gap. The Drude weight is given by 30-32 

1 ^ , 1 \- Kn^lJ^O;^)! 2 



1 9 2 -Eo(0) 

27V d<i>i 

= dH/d^, F^u 



(3) 

with J M = dH/d(j)^, F^y = d 2 H/d(j) ft d(f) u {n,u = x,y) 
and N being the number of sites. Here, E {(f>) and E n {4>) 
are the energy eigenvalues of the ground state |0; <f>) and 
of the n-th excited state |n; <f>), respectively. We use the 
second equation in cq. (3) to calculate the Drude weight, 
i.e., the second derivative of the ground-state energy with 
respect to the flux <p. We define the dimensionless Drude 
weight D = (D x + D y )l (47re 2 ). 33 

The charge gap is defined as a change in the chemical 
potential when electrons are added to or subtracted from 
the system, 

A c = li+-li-, (4) 



= -[E (N + 2,d>)-E (N,<f>)}, (5) 

M - = ±[E (N,<t>) -E (N-2,d>)]. (6) 

Here Eq(N b , <j>) is the ground-state energy with N e elec- 
trons, and the number of electrons is changed by two to 
conserve the z component of the total spin. 
We also calculate the double occupancy, 



docc = ^X^ n iT n a}<^ 



(7) 



which is believed to be related with the order parameter 
of the Mott transition. Here and henceforth (O)^ is un- 
derstood as the average value of the operator O in the 
lowest-energy state for the given <j>, (0)$ = (0; </>|O|0; <p). 

To study the magnetic properties, we calculate the spin 
structure factor, 



(8) 



We note that D, A c , d OCC} and S(q) are functions of 4>. 

3. Averaging over twisted boundary conditions 

In order to extract intrinsic properties of a bulk sys- 
tem, we need to extrapolate numerical results of finite- 
size clusters to the thermodynamic limit N — > oo. In 
practice, for small clusters that can be dealt with the ex- 
act diagonalization technique, the numerical results have 
large finite-size effects and strong dependence on the BC 
or 4>. Understanding this ^-dependence is an essential 
step towards systematic analysis of finite-size effects and 
extrapolation to the thermodynamic limit. 34 

To illustrate the problem we are facing, we first show 
^-dependence of the lowest energies of k = (0, 0) and 
(it, 7r) at U — 4 and 8 for t' = 0.5 on the 16-site clus- 
ter in Fig. 3. We clearly see that at U = 4 [Fig. 3(a)], 
the ground-state momentum changes from k = (0, 0) to 
(-7T, 7r) and vice versa as <fi varies. This implies that the 
states with different momenta become degenerate ground 
states in the thermodynamic limit, which is a typical fea- 
ture of a Fermi liquid (metal). On the other hand, as 
shown in Fig. 3(b), there is no momentum change in the 
ground state at U = 8, indicating the insulating behav- 
ior in the thermodynamic limit where a ground state is 
singled out with opening of a gap. 35 A similar behavior 
is seen already in the smallest 8-site cluster. 

Typically, there are three different ways to treat the 
BCs in the finite-size analysis: (i) fixing the BC, 36 i.e., fix- 
ing at a priori chosen value, for example, at (f) = (0, 0); 
(ii) optimizing <p for each cluster size and parameter set 
to minimize the ground-state energy, 37 ' 38 and; (iii) aver- 
aging over all the BCs computed physical quantities. 39,40 
The results obtained using these methods will, of course, 
converge to the same result in the thermodynamic limit, 
but can have quite different system-size dependences. 

The first approach (i) has ambiguity in choosing an ap- 
propriate BC for analyzing the size dependence, as finite- 
size effects appear very differently with different BCs. For 
example, when U is changed in the 16-site cluster, a level 
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Fig. 3. Boundary condition dependence of the lowest energies 
with k = (0, 0) and (tt, tt) at (a) U = 4 and (b) U = 8 for t' = 0.5 
on 16-site cluster. The inset in (a) shows the paramctrization of 

4>. 



crossing (or first-order transition) occurs for <p = (0,0), 
while no transition occurs for <fi = (n, 0). Thus, there is 
no obvious choice for a fixed and optimum BC, and we 
do not take this approach. 

The second approach (ii) also has a problem. For ex- 
ample, at U — the ground-state momentum changes at 
a finite value of t' that depends on the system size; see 
Fig. 4(a) below. Thus, it is difficult to obtain a smooth, 
systematic system-size dependence with this method. 

Compared with these two, the approach (iii), averaging 
over the BCs, does not have such difficulties and shows a 
well-behaved system-size dependence as we demonstrate 
below. In this approach, average is taken over the BCs 
for each cluster size, 39,40 



(O)bc = ± Jd<j> <0>«, 



(9) 



where A = (27r) 2 |bi x f»2 1 , and O is a physical quantity 
(such as the Drude weight D) introduced in the previous 
section. In the numerical calculations, the integral over 
the Brillouin zone is approximated by the summation 
over grid points. Error bars are estimated as a/ J 77$, 
where a is the standard deviation and N<p is the number 
of grid points. We take sufficiently large for 8- and 
10-site clusters. For 16-site cluster, we choose N$ — 32 ~ 
128 to obtain the precisions required. For 18-site cluster, 
we take N$ = 16. 

Let us demonstrate the efficiency of the averaging pro- 



cedure in comparison with the method (ii). Figure 4 
shows t' dependence of the Drude weight at U = ob- 
tained by the two methods (ii) and (iii). In Fig. 4(a) [the 
method (ii)], the Drude weight of each cluster size shows 
a sudden change (kink) in the t' dependence, which cor- 
responds to the change of the ground-state momentum. 
The kink appears at different values of t' for different sys- 
tem sizes, which complicates analysis of the system-size 
dependence in the method (ii). Furthermore, the results 
of finite N show large deviations of about 10% from the 
Drude weight in the thermodynamic limit. By contrast, 
Fig. 4(b) shows that the Drude weight averaged over the 
BCs [the method (iii)] at finite N is a smooth function of 
t' and is very close to that in the thermodynamic limit; 
The differences between the data points of N = 16 and 
N = oo are less than 1%. 

This good convergence, seen already at small N, can 
be understood as follows. At U = the Drude weight 
can be written in the thermodynamic limit as 



D = 



dk 



g(k)Q(e F - e k ), 



(10) 



where 0(x) is the Heaviside step function, cf the Fermi 
energy, and 

g(k) = 2t(cosk x + cosk y ) + At' cos(k x — k y ). (11) 

The single-particle energy is given by 

gfc = — 2£(cos£; a ; + cosfc y ) — 2t' cos{k x — k y ). (12) 

For a finite-size cluster the averaged Drude weight is cal- 
culated as 

{D) * C = 1 t\j ^ d ^H'9{k + 4>), (13) 

where the primed sum denotes summation over N lowest- 
energy single-particle states of energy e fc+( ^. We note that 
the number of single-particle states inside the Fermi sur- 
face determined by = £f is not always equal to N in 
finite-size clusters. Up to this finite-size effect eq. (13) 
is equal to eq. (10). This indicates that, in the non- 
interacting case, the thermodynamic limit can be reached 
by averaging over the BCs even in small clusters. 41 Note 
that the method (iii) can distinguish between a band 
metal and a band insulator. 

Next, Fig. 5 compares U dependences of the Drude 
weight at t' = obtained with the two methods (ii) and 
(iii) . This t' = case corresponds to the Hubbard model 
on the square lattice, in which the metal-insulator tran- 
sition is known to take place at U c — because of the 
perfect nesting. As shown in Fig. 5, however, the Drude 
weight computed with BCs minimizing the energy [the 
method (ii)] remains finite even when U t. A similar 
behavior was observed in the one-dimensional case. 34 In 
contrast to this, the Drude weight averaged over all the 
BCs [the method (iii)] immediately vanishes for U > 0, 
reproducing the expected metal-insulator transition, ex- 
cept for U « i where the averaged Drude weight will 
converge slowly to zero with N$ — > oo. 

The reason why the expected behavior is reproduced 
even with the small size clusters is understood as fol- 
lows. Let us consider the L x L clusters for simplicity. 
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Fig. 4. t' dependence of the Drude weight at U = obtained by 

(a) searching BC which minimizes the ground-state energy and 

(b) averaging over BC. 



At U = the ground-state momentum is k = (0, 0) for 
any </> because of the perfect nesting. In the </>i-<f>2 plane 
the ground-state energy has cusp singularities along the 
lines, 4>i — <p2 and 4>\ = —02, where the ground states 
are degenerate. As we discussed above, the Drude weight 
averaged over the BCs reproduces well the one in the 
N — > oo limit, when the cusp singularities are avoided in 
the average summation. A small but finite U lifts the de- 
generacy, and the cusp lines disappear. Then the ground- 
state energy becomes an analytic function of </>, and the 
average Drude weight vanishes by definition, for instance, 
(Ae)bc ot J dcf> d 2 E /d(j)1 = 0. This leads to the metal- 
insulator transition at U <C t expected for this case with 
the perfect nesting. 

From the results shown above, we conclude that the 
averaging over the twisted BCs is most efficient among 
the three methods to calculate D. We have applied this 
method to calculate the double occupancy and the spin 
structure factor, and confirmed its efficiency for these 
quantities as well. In the following section, we present 
numerical results of the Drude weight, the double occu- 
pancy and the spin structure factor, which are obtained 
with this method. (We will omit the brackets (• • • }bo) 

It turned out, however, that the method (iii) is not 
suitable for calculating the charge gap. To illustrate its 
drawback, let us consider a band insulator at U = 0. In 
this case, the charge gap is nothing but a band gap. If wc 
calculate the charge gap (4) from the chemical potentials 
/U± averaged over BC, the gap is overestimated, since the 
averaged chemical potentials [i + and /i~ correspond to 
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Fig. 5. U dependence of the Drude weight at t' = obtained by 
two different methods as in Fig. 4. 



the average energy at the bottom of a conduction band 
and that at the top of a valence band, respectively. To 
reproduce the band gap in the thermodynamic limit, we 
should instead search the minimum energy of the con- 
duction band and the maximum energy of the valence 
band. From this consideration, we use the definition 



A c = max { M + lin - Mmax , °} 



(14) 



where 



+ . E (N + 2, ( / ) )-E (N, 0) 
Mmin = mm 2 ' ( ) 



max ■ 



E (N,<j>)-E (N-2,<p) 



(16) 



With this definition, finite-size effects are completely re- 
moved for the band gap. Furthermore we obtain A c = 
at the noninteracting case U = as expected for a metal- 
lic state. These simple arguments suggest that the defi- 
nition of eq. (14) should be most effective for extracting 
the charge gap in the thermodynamic limit. Wc employ 
this definition in §4.2. 

4. Results 

In this section, we present our numerical results for D, 
A c , d occ and S(q) computed for a wide parameter range 
of t' and U with the method described in the previous 
sections. 

4-1 Drude weight 

Figure 6 shows the Drude weight at t' = 0.2, 0.5 and 
0.8. For each cluster size, the Drude weight decreases 
with increasing U and eventually vanishes. We can thus 
define the critical value Uj? above which the Drude 
weight vanishes. The critical value is expected to rep- 
resent the metal-insulator transition point. At t' = 0.2 
and 0.5, the system-size dependence of Uj? is small, and 
it is reasonable to estimate L/ C D = 3 ~ 4 for t' = 0.2 
and = 5^7 for t' = 0.5, respectively. However, 
at t' = 0.8 shows a substantial system-size dependence. 
Considering the fact that the Drude weight at U = 6 is 
almost independent of the system size, we can roughly 
estimate = 6^8. 

Let us remark on the discontinuity of the Drude weight 
at U — seen in the numerical results for small clusters. 
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Fig. 6. Drude weight at (a) t' = 0.2, (b) 0.5 and (c) 0.8 for 8-, 
10- and 16-site clusters. At t' = 0.5, data for 18-site cluster are 
also plotted. 

This discontinuity is an artifact of small system size and 
can be traced back to the cusp singularity at t' = 
mentioned in §3. In small clusters the number of discrete 
momenta is small, and the region represented by k = 
(0, 0) is relatively large. Hence, the effect of the collapse 
of the cusp lines becomes pronounced and leads to the 
discontinuity of the Drude weight even at finite t' . In 
fact, the jump in the Drude weight is larger for smaller 
clusters and for smaller t' as shown in Fig. 6. 

4-2 Charge gap 

Figure 7 shows numerical results for the charge gap 
A c defined in eq. (14). A finite t' lifts the perfect nest- 
ing, and the charge gap opens up when U is larger than 
some finite value. We thus define the second critical value 
above which we have a nonvanishing charge gap. Al- 
though it is difficult to see systematic size dependence, 



Fig. 7. Charge gap at (a) t' = 0.2, (b) 0.5 and (c) 0.8. 

we can roughly estimate the critical value as = 2^3, 
3.5 ~ 4.5 and 4 - 6 for t' = 0.2, 0.5 and 0.8, respectively. 

Note that the estimated values of are always 
slightly smaller than E/jP obtained from the Drude 
weight. 

4-3 Double occupancy 

Figure 8 shows the double occupancy d occ at t' = 0.2, 
0.5 and 0.8. We find that the size dependence of the dou- 
ble occupancy changes at some finite U. For example, at 
t' = 0.5, the double occupancy increases with increasing 
N for U < 5 while it decreases for U > 5. Accordingly, in 
larger clusters (e.g., N = 16), the U dependence has an 
inflection point. We thus define the third critical value 
lit? from the change in the size dependence. We estimate 
~ 4, 5 and 7 at t' = 0.2, 0.5 and 0.8, respectively. 

The double occupancy is considered to be closely re- 
lated with the order parameter of the Mott transition, 6 
and has been calculated to determine the transition point 




Fig. 8. Double occupancy at (a) t' = 0.2, (b) 0.5 and (c) 0.8. 

in the previous theoretical studies. We note that the val- 
ues of agree well with the critical values obtained 
from the Drude weight in §4.1. The relation among [/~ , 
C/jP and will be discussed in §5 in connection with 
the phase diagram. 

4-4 Spin structure factor 

In order to find possible magnetic phase transitions, we 
have calculated the spin structure factor, S(q), defined in 
eq. (8). For the parameters that we studied, S(q) always 
shows a peak at Q = (tt,tt), except when t' = 0.8 and 
U < 5 for the 16-site cluster. Thus we here analyze the 
staggered component, S(Q), only. 

Figure 9 shows the size dependence of S(Q)/N at t' = 
0.2, 0.5 and 0.8. Here we have also plotted the results 
for the 12-site clusters. Although there are two 12-site 
clusters with different geometries as shown in Fig. 2, the 
two results agree well, except at t' — 0.8. For t' = 0.2 
and 0.5, systematic size scaling is observed, whose scaling 



Fig. 9. Spin structure factor, S(Q) with Q = (tt,tt), at (a) tl = 
0.2, (b) 0.5 and (c) 0.8. The lines are fits by a linear function of 
7V-1/2. 

form is consistent with that expected in the staggered 
ordered phase, 

VS{Q)/N - moo oc TV-i/2 + . . . ) ( 17) 

where m^o — limjv— >oo y/S(Q)/N. From the value of U 
where moo becomes positive, we estimate the antifcrro- 
magnetic transition point to be at Uaf = 3.5 ~ 4.5 and 
5 ~ 6 for t' = 0.2 and 0.5, respectively. At t' = 0.8, 
however, we cannot find any systematic size dependence, 
suggesting that the antiferromagnetic state is not stable 
in this range of U. 

We note that the size dependences of S(Q) at t' — 0.5 
are consistent with the results obtained with the PIRG 
method for larger clusters N = 36 ~ 196. 29 This again 
confirms that the averaging over the twisted BCs cap- 
tures the asymptotic behavior in the limit of N — > oo 
from the calculations for small clusters. 
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5. Discussion 

In this section we discuss the numerical results pre- 
sented in §4 and deduce the ground-state phase diagram. 

5.1 Phase diagram 

In the previous section we have shown our numerical 
results of the Drude weight, the charge gap, the dou- 
ble occupancy and the spin structure factor, which were 
obtained from exact diagonalization of finite-size clus- 
ters up to 18 sites. Among them the Drude weight and 
the charge gap are the quantities directly signalling the 
metal-insulator transition. The critical values of U es- 
timated from these two quantities, however, turned out 
to be slightly different, < Uj? . On the other hand, 
U~ estimated from the double occupancy was about the 
same as l7 c d , i.e., Uj? ~ U^. Naively, these results would 
imply that, in the parameter region < U < U®, the 
charge gap opens but the Drude weight remains finite. 
It is, however, difficult to conceive a state with both a 
finite charge gap and a finite Drude weight, and we there- 
fore expect that, in the thermodynamic limit, the Drude 
weight should vanish whenever the charge gap opens. As 
discussed in the following, we speculate that the seem- 
ingly conflicting result < Uj? is due to non-trivial 
size dependence of D which cannot be captured by our 
calculations up to 18 sites. 

Here we propose two scenarios for the way in which 
our finite-size results approach the thermodynamic limit. 
The first scenario is obvious one: as the system size in- 
creases, the difference between Ujj" and Uj? is reduced 
and eventually vanishes, and we are left with one metal- 
insulator transition besides the magnetic transition at a 
larger value U = Uxe.\xi the other scenario there remain 
two phase boundaries corresponding to and Uj? even 
in the thermodynamic limit, and system-size dependence 
appears mainly in the Drude weight D (in the interval 
Uj^ < U < Uj?) which decreases to zero, while the dif- 
ference between U^ and Uj? is kept finite. 

First, let us consider the former scenario and discuss 
the size dependence of the Drude weight and the charge 
gap. Concerning the Drude weight, as we mentioned in 
§4.1, the data at t' — 0.8 clearly show a system-size de- 
pendence for U > 6, where the Drude weight becomes 
smaller in larger systems. Even in the case of t 1 = 0.5, 
if we look at the data in Fig. 6(b) carefully (by neglect- 
ing the results for the smallest 8-site cluster) , we can see 
smaller but similar system-size dependence for U > 4, 
although the variation of the data for different system 
sizes are within the error bars. We now show that these 
system-size dependences can be seen more clearly in the 
quantity which we denote by \D\, i.e., the average of the 
absolute value of the right-hand side of eq. (3) over the 
twisted BCs. The numerical results for \D\ are shown 
in Fig. 10. At t' = 0.5, \D\ converges on a single curve 
for U < 4 while it has systematic size dependence for 
U > 4. A similar size-dependent behavior is clearly seen 
for U > 5 at t' — 0.8. Although \D\ does not have a clear, 
direct physical meaning, these apparent system-size de- 
pendences of \D\ seen in Fig. 10 suggest that the Drude 
weight D retain a size dependence in the same range 
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Fig. 10. "Absolute value" of the Drude weight, \D\, at (a) t! = 0.5 
and (b) t' = 0.8. 



of U. In other words, the estimates of Uj? are consid- 
ered as upper bounds for the metal-insulator boundary 
monitored by the Drude weight. Note that the values 
of U where the size dependence of \D\ becomes evident 
roughly coincide with U^. 

As we have mentioned in §4.2, the charge gap does not 
show systematic size dependence up to N — 16. Inciden- 
tally, the PIRG study has shown that the charge gap of 
a 16-site cluster already gives a good estimate of that in 
the thermodynamic limit. 42 If we assume that this is also 
the case in the present results, then we may regard U^ 
obtained from up to 16-site clusters as a good estimate 
for the value in the thermodynamic limit. 

These considerations on the system-size dependences 
of the Drude weight and the charge gap led us to the first 
scenario in which, as N is increased beyond N = 18, Uj? 
decreases and finally coincides with Uj^ which changes 
little with N. In this case the Mott transition takes place 
at U^ . The insulating state with U larger than U^ does 
not show the antiferromagnetic order until U reaches 
Uaf which is substantially larger than U^. Although we 
do not have any further information on the magnetic 
properties of the insulating phase in the intermediate- 
coupling regime, one possibility is that the insulator is 
a featureless Mott insulating state without symmetry 
breaking, similar to the non-magnetic insulator found in 
the PIRG calculations. 29 Hence, in this scenario, we have 
three different phases in the plane of t' and U: a para- 
magnetic metal for U < U^, a non-magnetic insulating 
state for U^ < U < Uaf, and an antiferromagnetic in- 
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sulator for U > Uaf- 

Next, we explain the second scenario in which the 
phase diagram has two distinct phase boundaries cor- 
responding to Uj? and U^. In this case also, we as- 
sume that, in the thermodynamic limit, the system is 
insulating for U > because of the opening of the 
charge gap and because of possibly negligible size depen- 
dence of U^ as discussed above. One might then won- 
der how a finite Drude weight observed in the region 

< U < diminishes with increasing system size, 
while Uf? remains larger than U^. We speculate that 
this peculiar behavior of the Drude weight is a finger- 
print of some unconventional nature of the intermediate 
phase. Let us suppose, for example, that we have a bunch 
of low-energy levels which are separated by a gap from 
higher-energy modes. These low-energy levels are distin- 
guished by their total momenta k. Let us further assume 
that, as the system size increases, the band widths of the 
low-energy modes get narrower, while their dispersion 
curves remain to overlap as in Fig. 3(a). In this case we 
should observe finite Drude weight which is decreasing 
with increasing system size. When < U < , we 
thus expect to have an insulating phase with many low- 
energy modes, which is reminiscent of the non-magnetic 
Mott insulator found by the PIRG study. 29 Since the 
twisted BCs are not directly coupled to spin degrees of 
freedom, these low-energy modes involve some charge dy- 
namics and would not be obtained by spin-only models 
like the Heisenberg model with multiple-spin exchange 
interactions. The remaining problem is what happens at 
the phase boundary [/ C D as U is increased. A possibility 
is that there occurs a phase transition to another insu- 
lating state with some conventional symmetry breaking, 
e.g., an insulator with a complicated pattern of magnetic 
order or with dimer order. 43 ' 44 In such an insulating 
state, the ground state may consist of a single disper- 
sion as depicted in Fig. 3(b), and the Drude weight be- 
comes zero even in the small size clusters. Hence, in this 
second scenario, we have four different phases: a para- 
magnetic metal for U < U^, a non-magnetic insulating 
state for < U < U®, another insulating state for 

< U < Uaf, and an antiferromagnetic insulator for 
U > Uaf- 

We summarize the possible phase diagram in Fig. 11 
according to the two scenarios explained above. The grey 
curve represents the antiferromagnetic transition at U = 
Uaf- The dash-dotted curve corresponds to the opening 
of charge gap at U = U^. The broken curve denotes 
a possible phase transition from an unconventional, non- 
magnetic insulator to another insulating state at U = U® 
suggested in the second scenario. 

With the numerical results at hand, we cannot con- 
clude which scenario is more plausible. To answer this 
question requires further studies, such as calculations on 
larger clusters, a through examination of magnetic sym- 
metry breaking beyond the simple antiferromagnetic or- 
der with Q = (it, it) and of other symmetry breakings 
(dimerization etc.). 
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Fig. 11. Ground-state phase diagram of the Hubbard model on 
the anisotropic triangular lattice in t'-U plane. AFI and PM 
denote an antiferromagnetic insulating phase and a paramag- 
netic metallic phase. The dash-dotted curve represents the metal- 
insulator transition determined by the opening of the charge gap. 
The broken curve denotes a possible phase boundary between dif- 
ferent insulating states. The curves are the guides for eyes. See 
the text for details. 

5.2 Comparison with other theoretical studies 

Here, we discuss the implications of our results in com- 
parison with other previous theories, focusing on the na- 
ture of the Mott transition. Comparing the phase dia- 
gram in Fig. 11 with that obtained by the large-scale 
PIRG calculations, 29 we find that the Mott transition 
line in our phase diagram is in good agreement with the 
one determined by the PIRG method. However, there is 
an important discrepancy in the evolution of the double 
occupancy as a function of U. 

In our results, as shown in Fig. 8, the double occupancy 
is a smooth decreasing function of U, showing no notable 
anomaly in the Mott transition at U^ . On the contrary, 
in the PIRG results, it shows a jump at the transition, 
indicating the first-order nature of this metal-insulator 
transition. Moreover, in the metallic state in the vicinity 
of the transition, the double occupancy is not a decreas- 
ing function of U; it is almost flat in a certain range of 
U and keeps a relatively large value of about 0.2. This 
flat regime as well as the jump is not reproduced in our 
calculation of the double occupancy. 

On the other hand, in the insulating side just above 
l7 c a , the double occupancy shows a good agreement be- 
tween the two results. Our results indicate an opening 
of the charge gap expected there. Furthermore, as we 
discussed in §5.1, a finite Drude weight showing a small 
system-size dependence can be ascribed to some uncon- 
ventional nature of the insulating state, as suggested in 
the PIRG study. 29 Hence we consider that our present 
scheme has a potential to describe the insulating side of 
the Mott transition. 

Thus the most peculiar feature is seen in the metallic 
side of the Mott transition. There the PIRG results indi- 
cate that carriers have coherent nature even at a moder- 
ately large value of U. In fact, the momentum distribu- 
tion shows a large jump at the Fermi surface without sig- 
nificant reduction of the renormalization factor. 29 If we 
regard the PIRG results as standard results to be repro- 
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duced, then the discrepancy between our results and the 
PIRG results implies that the coherent nature of a carrier 
is not captured by our calculations for small size clusters, 
even though strong correlations are taken into account 
properly in our exact-diagonalization study. Moreover, 
the peculiar behavior of the double occupancy is not 
reproduced by both the VMC study 22 and the cellular 
DMFT study. 23 Considering the fact that the system size 
used in the VMC study is much larger than the one in 
our study, we suspect that the discrepancy between the 
VMC study and the PIRG is indicating that the many- 
body effects relevant to the peculiar coherency are not 
treated properly in the VMC which employs approxi- 
mation of a single Slater determinant with a Gutzwillcr 
projection. The disagreement between the results of the 
cellular DMFT and the PIRG might be ascribed to the 
problem shared by our study, i.e., lack of sufficient infor- 
mation on the spatial fluctuations in larger scale. Further 
study is highly desired to clarify the peculiar nature of 
the metallic state, and this will help us better understand 
the Mott transition. 

6. Summary 

We have studied the ground-state properties of the 
two-dimensional Hubbard model on the anisotropic tri- 
angular lattice by the exact diagonalization method. In 
order to reduce the finite-size effects, we have demon- 
strated efficiency of averaging over the twisted boundary 
conditions for reducing finite-size effects. Using this tech- 
nique, we have calculated the Drude weight, the charge 
gap, the double occupancy and the spin structure factor, 
and determined the the phase boundaries of the metal- 
insulator transition and the antiferromagnetic transition. 
An unconventional aspect of the Mott transition of the 
present system is illuminated by comparing our results 
with other previous theoretical results. 
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